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Abstract 

Measurements of the cosmic-ray air-shower fluorescence at extreme energies require pre- 
cise knowledge of atmospheric conditions. The absolute calibration of the cosmic -ray en- 
ergy depends on the absorption of fluorescence light between its origin and point of its 
detection. To reconstruct basic atmospheric parameters we review a novel analysis method 
based on two- and multi-angle measurements performed by the scanning backscatter lidar 
system. Applied inversion methods, optical depth, absorption and backscatter coefficient, 
as well as other parameters that enter the lidar equation are discussed in connection to the 
attenuation of the light traveling from shower to fluorescence detector. 

Key words: backscatter lidar, inversion methods, two- and multi-angle reconstruction, 
atmospheric optical depth, cosmic showers, fluorescence detectors 
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1 Introduction 



Contemporary fluorescence experiments (Fly's Eye [1], HiRes [2], P. Auger [3]) 
studying cosmic rays with energies near 10 20 eV detect fluorescence light produced 
along the air-shower volume. As a cosmic -ray induced particle cascade develops 
in the atmosphere it dissipates much of its energy by exciting and ionizing air 
molecules. The cascade particles along the shower axis are limited to a narrow 
lateral distribution where excited nitrogen molecules are fluorescing in the near- 
UV spectral band. The efficiency of the process is rather small. However, the vast 
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number of emitting particles in high energy shower makes this source of radiation 
highly significant. Ultimately, this electromagnetic (EM) cascade dissipates much 
of the primary particle's energy. Fluorescence light is emitted isotropically with an 
intensity proportional to the number of charged particles in the shower. EM com- 
ponent and hence the total number of low energy EM particles is in turn fairly 
accurately proportional to the energy of the primary particle. Thus, the calorimetric 
measure of the total EM shower energy [4] is proportional to the integral of EM 
particle density iV em along the shower direction x, 



E em = K 



J N em (x) dx (1) 



with K ~ 2.2 MeVcm 2 /g, where x is measured in units of longitudinal air density 
(g/cm 2 ). E em is a lower bound for the energy of the primary cosmic ray. The lower 
portion of shower development is usually obscured by the ground so that EM cas- 
cade reaching below ground is included by fitting a functional form to the observed 
longitudinal profile and integrating the function past surface depth. The number 
of photons iVph reaching fluorescence detector (FD) is proportional to EM particle 
density N em (x) at the point of production x, so that in turn 

T(x) 

with R(x) being distance between shower point x and FD. Light originating within 
the shower is certainly affected by the absorption and scattering on molecules and 
aerosols in the atmosphere. The number of detected photons is thus reduced due to 
non-ideal atmospheric transmission T(x) < 1, where 



T(x) = exp 



— / a(r) dr 
o 



e~ T{x \ (3) 



with a(r) volume extinction coefficient along the line-of-sight, and t(x) the result- 
ing atmospheric optical depth (OD) to the shower point x. 

In this sense, the atmosphere can be treated as an elementary-particle detector. 
However, weather conditions change the atmospheric transmission properties dra- 
matically resulting in strongly time-dependent detection efficiency. Therefore, an 
absolute calibration system for fluorescence light absorption is an essential part of 
FD [5,6]. 

Eq. (3) is a basis for the fluorescence-detector energy calibration. Apart from the 
Pierre Auger Observatory, all existing fluorescence-based experiments have suf- 
fered from a lack of a sufficient atmosphere monitoring system. Weather conditions 
in a desert-like atmosphere were expected to be stable enough, so that the standard 
attenuation-length profile should suffice to reconstruct the total EM shower en- 
ergy, Eq. (1), with a controllable precision. However, this has turned out not to be 
the case, especially for rare events with energies above 10 19 eV. In addition, the 
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energy reconstruction is obscured by Cerenkov radiation in the lower part of the 
air-shower which can not be separated from fluorescence. For more than a half of 
the highest-energy events measured so far the atmosphere properties are not known 
well enough to accurately reconstruct the primary energy. To be able to provide 
adequate calibration, one has to measure attenuation at the time of the event in the 
whole region of the air-shower. 

There is also a systematic discrepancy when comparing the cosmic-ray spectra of 
fluorescence experiments with ground arrays. Their compatibility can be estab- 
lished only with a shift in energy and flux of one or the other. As in the case of flu- 
orescence detectors, ground arrays have their own problems with the energy deter- 
mination and are much more dependent on air- shower simulations. At present, it is 
not known whether the discrepancy is due to fluorescence-detector or ground-array 
method, or both. Therefore, it is of utmost importance to have in situ atmosphere 
monitoring system which is working coherently with a fluorescence detector. 

To lower primary cosmic ray energy uncertainties, the volume extinction coefficient 
a(r) thus has to be well estimated over almost whole detection volume of FD. In the 
case of the Pierre Auger Observatory, the detection volume corresponds to ground 
area of 3000 km 2 and height of ~ 15 km. 

The paper is organized as follows. The first two sections are devoted to introduc- 
tory material on lidar measurement technique and a description of our specific ex- 
perimental setup. Then the atmospheric model for simulation of lidar signals is 
presented and the signals are evaluated by two well established inversion meth- 
ods. Results of the inversions are compared to the input model and conclusions on 
their applicability are drawn. Next, improved approaches to FD calibration based 
on scanning lidar system are proposed and evaluated on real data obtained with our 
experimental setup. 



2 Lidar system 



One of the most suitable calibration setups for FD is the backscattering lidar system, 
where a short laser light pulse is transmitted from FD position in the direction of 
interest. With a mirror and a photomultiplier tube, backscattered light is collected 
and recorded as a function of time, i.e. as a function of backscatter distance. Note 
that light from the lidar source traverses both directions, so that in case of matching 
laser and fluorescence light wavelength, OD for lidar light sums to twice the OD for 
fluorescence. The lidar equation [11] describes the received laser power P{r) from 
distance r as a function of volume extinction coefficient ct(r) and backscattering 
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coefficient @(r), 

P(r) = P ^(r)4e- 2rW . (4) 

P is the transmitted laser power and A is an effective receiving area of the detector, 
proportional to the area of the mirror and proportional to an overlap between its 
field of view with the laser beam. t is laser pulse duration. As seen from Eq. (2), 
measurement precision of a and corresponding r directly influences the precision 
of primary particle energy estimation. 

Simple as it may look, the lidar equation (4) is nevertheless difficult to solve for 
two unknown variables, ct(r) and (3{r). All existing analysis algorithms (Klett [7], 
Fernald [8], and their respective variations) reviewed in one of the following sec- 
tions are based on an experimental setup with static beam direction. This leads to 
ambiguity in determination of a(r) and f3(r) which can not be resolved without 
additional assumptions about atmospheric properties. At the FD experimental sites 
the atmosphere can be assumed to be almost horizontally invariant. In this case, 
there is an additional constraint when comparing signals coming from different di- 
rections, which solves the lidar equation for a(r) and (3{r) unambiguously. Even 
at this point, it can be stated that the need for steerable (scanning) lidar setup is 
unavoidable for proper solution of lidar equation. 



3 Experimental setup 

The lidar system used for verification of the analysis method is based on the Contin- 
uum MiniLite-1 frequency tripled Nd:YaG laser, which is able to transmit up to 15 
shots per second, each with energy of 6mJ and 4 ns duration (1.2m). The emitted 
wavelength of 355 nm is in the 300 — 400 nm range of fluorescence spectrum. The 
receiver was constructed using 80 cm diameter parabolic mirror with focal length 
of 41 cm. The mirror is made of aluminum coated pyrex and protected with Si0 2 . 

The backscattered light is detected by a Hammamatsu R7400 photomultiplier with 
operating voltage up to 1000 V and gain of 10 5 to 10 6 . To suppress background, 
a broadband UG-1 filter with 60% trasmitance at 353 nm and FWHM of 50 nm is 
used. The distance between laser beam and the mirror center is fixed to 1 m, and 
the system is fully steerable with 0.1° angular resolution. 

The signal is digitized using a three-channel LICEL transient recorder TR40-160 
with 1 2 bit resolution at 40 MHz sampling rate with 1 6k trace length combined with 
250 MHz photon counting system. Maximum detection distance of the hardware is 
thus, with this sampling rate and trace length, set to 60 km. However, in real mea- 
surements, atmospheric features up to 30 km only are observed. LICEL is operated 
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Figure 1. Schematic view of the lidar system. A mirror of 80 cm diameter and a UV-laser 
head are mounted on the steerable mechanism. The LICEL TR40-160 receives the trigger 
from the laser and the signal from Hammamatsu R7400 phototube. The Linux-PC controls 
the LICEL digitizer through PCI-DIO-32HS Digital Input/Output card. The steering motors 
are controlled through RS-232 port. Zenith angle is denoted by 0. 

using a PC-Linux system through a National Instruments digital input-output card 
(PCI-DIO-32HS) with Comedi drivers [9] and a ROOT interface [10]. 



4 Lidar simulation with specific atmospheric model 



In a low opacity atmosphere the attenuation and the backscattering coefficient can 
be written as a sum of contributions from two independent components, 

a(h) = a m (h) + a a (h), (5a) 
P(h) = P m (180> m (/i) + P a (180> a (». (5b) 

where a m and a a correspond to molecular and aerosol attenuation, respectively. 
The aerosol phase function P a (180°) for backscattering has, apart from the wave- 
length, also a strong dependence on the optical and geometrical properties of the 
aerosol particles. Nevertheless, at wavelength of 355 nm, values in the range 0.025 
and up to 0.05 sr _1 can be assumed [11] for aerosol phase function P a (180°). The 
angular dependence of molecular phase function is defined by the Rayleigh scatter- 
ing theory, where P m (180°) = 3/8tt sr" 1 . 

For simulation purposes, the elevation dependence of the extinction coefficients is 
modelled as following, 

a m (h) 
a a (h) 

where L m and L a are the molecular and aerosol attenuation lengths at ground level, 
and and /i° are the molecular and aerosol scale height, respectively. An addi- 
tional mixing height h x is set up for aerosols, assuming uniform concentration near 
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Figure 2. Extinction-backscatter plot (a(3 diagram) for the model atmosphere in Eq. (6). 

the ground level and continuous transition into exponential vanishing for h > h x . 
The following values of the parameters are used: L m = 15 km, = 17.5km, 
L a = 2 km, h x = 0.8 km, and h° = 1.4km. 

The atmospheric model, Eq. (6), serves as a testing ground for two widely used 
reconstruction methods presented in the next section. A comparison with recon- 
struction of the real atmosphere yields insight into the common problems of the li- 
dar field. The Poissonian statistics of photon counting and multiplying, background 
noise, and effects of digitalization have been taken into account in the generation 
of the simulated lidar signals and under inspection match those observed in the real 
lidar power returns. 

Model in Eq. (6) is a valid approximation to the atmospheric conditions found in 
real experiments. Although the vertical variation of aerosol and molecular densi- 
ties is quite simple, the model still produces a non-trivial relation between total 
attenuation a and total backscattering coefficient (3. Therefore, the dependence of 
f3 as a function of a, shown in Fig. 2, cannot be well approximated by some simple 
functional form. 




5 Reconstruction of a ID atmosphere 



Concentrating on a single shot lidar measurement, the optical properties obviously 
have to be reconstructed in a ID subspace of the atmosphere. Rewriting the lidar 
equation (4), 

P( r ) = jB ^l e - 2 ^) (7) 
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where the effective aperture of the system is gathered in the constant B, an auxiliary 
S-function can be introduced, 

S (r) = In ^f^- 2 = In [0(r)/A>] - 2r(r; r ). (8) 

Note that r(r; r ) = J r r a(r') dr' corresponds to atmospheric OD between r and 
r. 



5.1 Klett inversion 



Apart from the experimentally measured lidar power return P(r), in Eq. (7) there 
are two unknown quantities, j3 and a (or equivalently r), preventing the unique 
solution of the lidar equation. Nevertheless, a simple, and sometimes physically 
meaningful, assumption of proportionality between backscattering and extinction, 

P(r) oc a k (r), (9) 

allows for the transformation of the integral Eq. (8) to the corresponding Bernoulli's 
differential equation with an existing analytical solution. Direct application of the 
solution (forward inversion) is numerically unstable, in some cases singular, and 
highly sensitive to the signal noise [7,12]. Klett's reformulation [7] of the solution 
(backward inversion) avoids these problems. The lidar backward inversion algo- 
rithm proceeds from the far point of the measured signal r f to the near end, 

e S(r)/k 

a(r; ai) = e*/*/a f +f /; f eS(*-')/* dr'' (10) 

where Sf = S(r f ), and ctf = a(r f ) is an estimate for the attenuation at the far 
end of the data set. The reconstructed attenuation a(r; a f ) is still a one-parameter 
function of the unknown boundary attenuation value « f , so that independent mea- 
surement, or suitable approximation, is needed at the reference distance r f . OD can 
be expressed directly from Eq. (10), 



k 

r(r; r , a f ) = - In 



fce 5f / fc +2a f / r r o f e 5 ( r ')/ fe dr' 



ke s ^ k +2a f j; { e s ( r ')/ k dr' 



(ID 



Klett's inversion method depends rather strongly on the assumed power law pro- 
portionality in Eq. (9). In Fig. 3, a failure of this approximation is demonstrated 
for the specific atmospheric model used for our simulations. The local value of the 
exponent, 



a dp 
= ^Ir" 



da 
dr 



(12) 
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Figure 3. Effective power k in Eq. (9) as obtained from the model atmosphere in Eqs. (5) 
and (6). Note that the discontinuity at h x = 1.4 km arises due to aerosol part of the model, 
up to where aerosol concentration is kept constant. 



Figure 4. Reconstructed attenuation a(h) from Klett's inversion of the simulated vertical 
shot data as obtained by different boundary values ctf. Solutions with 0.5, 1, and 2 times 
the correct Qf are plotted with dots. The actual model a profile is drawn with solid line. 
Assuming range-independent (constant) Klett's k, the best agreement between the recon- 
structed and actual profile is achieved for k 0.5, therefore this value is used for all tree 
plots. 

is shown to possess substantial range dependence. The main reason for failure of the 
power law proportionality stems from the inequality of the molecular and aerosol 
phase functions, P m (180°) and P a (180°), rendering the a and (3 relationship de- 
pendent on the particular magnitude of both quantities, and consequently range 
dependence (see Fig. 2). Therefore, the best value of k must be chosen using some 
ad hoc criterion. 

Analyzing results in Fig. 4, presenting Klett inversion of simulated lidar signals, it 
seems that the closest reconstruction of the model profile is achieved with k ps 0.5. 
From Fig. 3, showing local exponent k obtained with use of Eq. (12), it can be seen 
that k 0.5 is observed only in small interval around 4 km whereas at other places 
it is substantially larger. For r > 8 km, dominated by molecular scattering it, slowly 
approaches the value of 1, most commonly adopted in the literature. Nevertheless, 
as can be seen in Fig. 5, reconstruction of OD with k — 1 totally fails to reproduce 
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Figure 5. Reconstructed optical depth r(h; ho) as obtained from Klett's inversion in Fig. 4 
(with 0.5, 1, 2 times the correct boundary ccf). Upper three curves, with different ccf, are 
obtained using fc = 0.5 and lower three with the most frequently used k = 1. The actual r 
profile from our model is drawn with solid line and is seen to be well approximated by the 
middle curve with k = 0.5. 

the correct answer. Surprisingly, in case of this specific atmospheric model the most 
authentic result is obtained with k 0.5. 

Another drawback of the Klett's method is estimation of the extinction ctf at the far 
end of the lidar return. In the case that r t corresponds to a highly elevated point, ap- 
proximation a f = a m (r f ), i.e. the extinction at that point is dominated by the molec- 
ular scattering, yields quite reasonable results [13] with qualitative convergence to 
the correct a-profile. In general, for optically dense atmosphere (e.g. presence of 
moderate haze) convergence of the Klett's method is far more rapid as in clear, opti- 
cally thin case. However, sites for FD are usually chosen at locations with clear and 
cloudless atmosphere. For horizontal lidar measurements (zenith angle = 90°) in 
a horizontally invariant atmosphere, a t can be estimated as the one that minimizes 
extinction deviations from a constant value [13,14], i.e. minimizes the functional 



5.2 Fernald inversion 



Since the concentration of the molecules depends solely on the thermodynamic 
parameters of the atmosphere, the Rayleigh scattering on molecules is modeled 
separately on a basis of the meteorological data. ct m (r) acquired in that way is 
inserted in Eq. (5). With an estimate for the molecular and aerosol backscattering 
phase fraction, F = P m (180°)/P a (180°), and modified S-function 



J>(r')-a f ] 



2 dr'. 




(13) 
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Figure 6. Fernald inversion of simulated lidar signal. The correct result is drawn in solid 
line. The three data sets are inversions with a a (ri) = and ±a m (Vf)/2. The phase fraction 
F is kept equal to the value used for generation of simulated lidar returns. 

the lidar equation can be solved for aerosol part a a (r) following the same steps as 
in Klett's version, 



a a (r) = -Fa m (r) + 



(14) 



with ctf = Fa m (rf) + a a (r f ) and Sf = S(rf) = S(r f ). In the same way OD is 
expressed as 



r(r;r ,a f ) 



iln 

2 



e*+2af L rm e^ r ') dr' 

1 

e^+2a f /; m e^') dr' 
+ (1 - F) f r a m (r') dr'. 



+ 



(15) 



Note that the Fernald procedure relies on three independently supplied parameters: 
(i) an accurate estimate of molecular part of the scattering a m (r) along the whole 
range of interest, (ii) total extinction at the far end a f , and (iii) proper approximation 
for phase fraction F. As predicted by the Mie theory, it is quite difficult to obtain 
reasonable values for the latter. As for <S t , conclusions are similar to those of Klett's 

Ctf. 



In Fig. 6 Fernald's inversion of simulated lidar return is shown for different input 
values of a a (rf) that enter total extinction a t . For upward pointing lidar measure- 
ments vanishing aerosol concentration can be assumed at the far end of atmosphere, 
i.e. a a (r f ) = 0. To test the sensitivity of the reconstructed OD on this assump- 
tion, data sets with a a (r f ) = ±a m (r f )/2 and therefore <S f = (F ± l/2)a m (r f ), are 
also plotted. P a (180°) = 0.025 sr -1 is used in the expression for phase fraction F. 
Comparing to the Klett's method, which does not separate aerosol and molecular 
scattering, it is not surprising that the variation of Fernald's results on boundary pa- 
rameters is somewhat weaker. Pinning the molecular part of scattering undoubtedly 
stabilizes OD profiles obtained. Nevertheless, Fernald's inversion still relies heav- 
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ily on additional external parameters that are usually difficult, if not impossible, to 
measure. 



6 Horizontally invariant atmosphere 

Fluorescence detectors for cosmic showers are usually placed at locations with spe- 
cific atmospheric conditions. In case of the Pierre Auger Observatory, the FD cam- 
eras are covering the lower part of the atmosphere over an almost perfect 3000 km 2 
plane 1500 m above the sea level with remarkable fraction of cloudless days. Due 
to the high elevation and dry inland climate, an optically thin atmosphere is ex- 
pected. But, as noted before, in this case convergence of Klett's method is slower 
and can lead to erroneous estimates of OD. Based on that, and other peculiar prob- 
lems of the well established lidar inversion methods, a new approach with fewer 
a priori or hard-to-estimate input parameters is needed. Since the lidar equation 
is not uniquely solvable, a minimal set of assumptions needed for inversion has to 
be reconsidered. For a typical FD site it is quite reasonable to assume weak hor- 
izontal variation of the atmospheric optical properties. That is even more true for 
the huge plane mentioned above, with hardly any changes in elevation and vege- 
tation coverage. Since the FD is exclusively operating at night, only atmospheric 
conditions at that time have to be considered. The mean night wind speeds do not 
exceed 12km/h [15], so that particularly thin layer of aerosols close to the ground 
is expected. At night, it is also expected that there will be a low probability for 
formation of convective types of atmospheric instabilities. 

6.1 Two-angle reconstruction 

Under the moderate assumptions presented above, optical parameters of atmo- 
sphere that enter the lidar equation (7) can be assumed to possess only vertical 
variations, while being uniform and invariant in the horizontal plane. 

Thus, it makes sense to rewrite the range dependent S'-function in Eq. (8) in terms 
of height h and geometric factor £ = 1/ cos = sec 0, when lidar shots with zenith 
angle are considered. The S'-function becomes 



with "vertical" OD r(h; h ) = f£ Q a(ti) dti and (3 = /3(h ). After measuring two 
S-functions at different zenith angles £i = 1/ cos 0i and £2 = 1/ cos 2 and height 
h, Eq. (16) can be solved for the vertical OD, 



S(h,£) = hi\p(h)/Po]-2ZT(h; h ) 



(16) 





(17) 



2 



6-6 
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and the backscatter coefficient ratio, 



A) 



= exp 



6-6 



(18) 



Both quantities are directly proportional to the difference of two ^-functions at the 
same height and different angles. Therefore, choosing a small separation between 
zenith angles, 6 = £ an( l 6 — £+ d£, a differential form of Eq. (17) can be written, 



r{h) = 



ldS_ 
'2d£ 



(19) 



Equivalently, the differential form of Eq. (18) can be obtained, 



Po 



exp 



(20) 



Note that the OD is in that way determined up to the additive constant, and the 
backscatter coefficient up to the multiplicative factor. Nevertheless, both values 
should satisfy S(h ) = and r(h ) = 0. 

Taking into account the Poissonian statistics of collected photons, and neglecting 
all other sources of measurement uncertainties, a relative error of the obtained OD 
at some height depends on the lidar system parameters, 



h/ho 



16 - 6 



e 2£ir _|_ e 2^ 2 T ) 



(21) 



' 2t^N [3 
as well as relative error of backscatter coefficient 



P 



h/h y/&e&T+&&*T 



N P 



16 - 6 



(22) 



where N is number of detected photons in the time interval corresponding to the 
power return from height h , and f3 = P/Po- 

In Fig. 7, an example of ^-functions and their zenith angle variation is presented. 
All results are obtained from real lidar measurements performed during few Novem- 
ber nights in a typical urban atmosphere (GPS location: 46°04'35" N, 014°29'05" 
E, 312 m above sea level). For fixed primary azimuth angle <j>\ = 0° and three 
selected secondary angles 2 = 38°, 42°, and 47° results for OD (Fig. 8), backscat- 
ter coefficient (Fig. 9), and aP diagram (Fig. 10) are obtained from corresponding 
S-f unctions in Fig. 7. Due to presence of a thin layer of optically thick haze at 
h ~ 3 km, a drastic change in both OD and backscattering at that height is ob- 
served. Since OD is well determined only up to an additive constant, note that the 
variation of results for different 2 is easily produced by the inadequate determi- 
nation of So, in other terms, by variation of atmospheric optical properties at ho. 
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Figure 7. S-function at few angles, 4> = 0°, i.e. £ = 1 (upper data set in black), 
42°, and 47° (£ = 1.27, 1.35, and 1.47), in shades of gray (lower three data sets). 
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Figure 8. Reconstructed optical depth (OD) r from three pairs of S-functions. In all pairs, 
Si corresponds to the S-function with cj) = 0° (£ = 1) and 52 to the S'-functions with 
(j) = 38°, 42°, and 47° (£ = 1.27, 1.35, and 1.47), respectively from bottom to top. 
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Figure 9. Reconstructed backscatter coefficient /3(h) /(3(ho) from three pairs of 
S'-functions. In all pairs, Si corresponds to the S-function with cj) = 0° (£ = 1) and S2 
to the S-functions with = 38°, 42°, and 47° (£ = 1.27, 1.35, and 1.47), respectively 
from bottom to top. 

Compatible with a scale height of ~ 18 km, the variation of backscattering in Fig. 9 
is slower as found in our model, generating a gradual but still comparable aj3 dia- 
gram in Fig. 10. 
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Figure 10. a[3 diagram (extinction-to-backscatter plot) for the pair where S\ is taken at 
4> = 0° and S2 at <p = 38°. Note that in Fig. 2 the same diagram is plotted also for model 
atmosphere. 
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Figure 11. Logarithmic plot of the relative deviation in OD, <j t /t (upper panel), and 
backscattering coefficient, <Jg//3 (lower pannel), vs. the second shot angle </>2, in case that 
the first shot zenith angle is set to <\>\ = 0°. Values h/ho = 8, 7V = 4 ■ 10 6 , r = 1, and 
P/Po = 0.6 corresponding to the far point (h w 8 km) in Fig. 8 have been assumed for 
parameters in Eqs. (21) and (22) (solid line). Values h/h = 3, 7V = 4 ■ 10 6 , r = 0.4, and 
P/Po = 0-8 corresponding to the near point (h ps 2 km) are assumed for the dashed curves. 
Note that <\> = 60° corresponds to £ = 2. 

In Fig. 11, a logarithmic plot of the relative error in OD is presented for typical 
lidar system parameters. First, the angle is fixed to 0i = 0° while the second one, 
02, is varied from a vertical to an almost horizontal shot. It is hard to avoid the fact 
that minimum error is produced with evaluation of two quite considerably separated 
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lidar shots, </> 2 ~ 70°. Even at moderate elevations h this can amount to large spatial 
separations of the two points of lidar return, and thus the requirement of horizontal 
invariance easily broken. In the case of an atmosphere, that is slowly horizontally 
modulated, a more "local" approach to the OD problem is needed. 

6.2 Multi-angle reconstruction 

For the ideal atmosphere, with true horizontal invariance, the £ dependence of S- 
function is particularly simple, 

S(h,0 = ln[P(h)/(3 ]-2^r(h;h ), (23) 

with the backscatter coefficient m[/3//3o] as offset, and OD r as the slope of the re- 
sulting linear function in £. Therefore, the optical properties of the atmosphere can 
be alternatively obtained from the analysis of the ^-function behavior for scanning 
lidar measurements. Furthermore, disagreement of the measured S(£) profiles from 
the linear form is a suitable criterion for detection of deviations from the assumed 
horizontal invariance of the atmosphere. 

A generalization of the two-angle equations (19) and (20) to their differential coun- 
terparts strongly suggested this way of reconstruction of optical properties, the 
two-angle method being a mere two-point approximation of the linear function 
in Eq. (23). Taking into account quite substantial uncertainties in S(£) for single 
angle, the linear fit trough many data points seems to yield superior results and 
the reconstruction is no longer limited to two lidar shots, well- separated in angle. 
The preferred horizontal invariance is not required to take place across huge atmo- 
spheric volumes (as in case of 0i = 0° and 2 = 60° shots), but has to be met only 
in relatively small arc of interest where the continuous lidar scan is performed. 

In the opposite case, when slow variation of atmospheric properties in horizontal 
plane is allowed, Eq. (19) is similar enough to the renown ID "slope method", 
based on assumption of small variation of /3(r), or equivalently d/3/ dr 0. Bear 
in mind that in method presented here the variation of (3 with height can be of any 
magnitude, as long as there are only modest variations in the horizontal direction. 

Opposite to Fig. 7, in Fig. 12 ^-function profiles with respect to zenith £ are drawn 
for fixed heights, starting with h = 3.2 km and up to 7 km with 633 m step. Approx- 
imate linear behavior is observed in few arc intervals, with narrow bands of minute 
atmospheric shifts at £ = 1.15 and 1.38. Since these shifts in profiles disappear 
when lifting h from 3 km to 3.5 km, they are obviously due to the distortions of 
atmosphere in the latter interval, feature already observed in Fig. 8. 

In Fig. 13 the results of fitting and extraction of OD are similar to the ones in Fig. 8. 
Note that in both cases OD is obtained relative to the h point, so that the results 
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Figure 12. Dependence of S'-function on azimuth angle at various heights h = 3.2 (black), 
5.6 (gray), and 8 km (light gray), while ho = 3km. Note that £ = 1 corresponds to cf) = 0°, 
and £ = 1.5 to = 48°. 
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Figure 13. Optical depth r obtained by linear fits of angle dependence of S-f unctions in 
Fig. 12. 

may differ up to some additive constant. Therefore, comparing both figures, it is 
more accurate to concentrate on the same span of OD within the 3.5 km to 9 km 
interval. Nevertheless, the range of OD results with acceptable error bars is with 
multi-angle method increased up to 12 km. 

The relative error of OD in Fig. 14 is needed for correct estimation of shower energy 
uncertainty. It is kept below 6% even for the OD from the far points of the range, 
and below 3% for modest values of OD. Fig. 15, with values for (3/{3 should be 
compared to Fig. 9. 



7 Conclusions 



Inversion attempts of simulated lidar returns for atmosphere, modeled by Eqs. (6), 
show numerous drawbacks of established numerical methods. For instance, Klett's 
and Fernald's method of section 5.1 and 5.2 do not satisfy the specific requirements 
of FD calibration. While they may be useful for qualitative reconstruction of atmo- 
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Figure 14. Dependence of relative error in optical depth on depth itself, a T /r. Data points 
and uncertainties are from Fig. 13. 
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Figure 15. Relative backscattering coefficient (3(h) /(3q from S'-functions in Fig. 12. 

spheric properties (spatial haze/cloud distribution, cloud base etc.), they are not ap- 
plicable for absolute assessment of atmospheric transmission properties. There are 
many reasons for this failure. One of them is certainly strong dependence of ob- 
tained inversions on presumed extinction/backscatter functional relation, Eq. (9), 
in case of Klett's method, and assumed spatial dependence of Rayleigh scattering 
on molecules in Fernald's case. Another issue is the extraordinarily difficult mea- 
surement of far- side extinction rate ctf, needed in Eq. (10), and phase fraction F, 
Eq. (13). We are therefore forced to find better solutions, even at the expense of 
adding scanning capabilities to an otherwise rigid lidar setup. 

In contrast to that, based on the sole assumption of a horizontally invariant (or at 
least horizontally slowly varying) atmosphere, the two-, and especially the multi- 
angle, method presented in section 6, while simple in structure, nevertheless pro- 
duce reliable quantitative answers with small uncertainties (e.g., see Figs. 8 and 
13) to FD calibration questions. As found by our investigation of first-run measure- 
ments, the relative error of OD for distances up to 12 km stay well below 6%. This 
number can be reduced even further by slow angular scanning and fast multiple- 
shot averaging of lidar returns. Nevertheless, in that case increased interaction be- 
tween FD and lidar laser source, especially FD blind time, has to be taken into ac- 
count. Furthermore, concerning the specific form of the atmospheric transmission 
entering Eq. (2), they offer suitable starting ground for development of methods 
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that can considerably reduce systematic errors of shower energy E em estimation 
with fluorescence detectors. 

In the case of strict horizontal invariance, both methods deliver exact solutions of 
the lidar Eq. (4) with accuracy of the results limited only by the quality of the 
measurement. In that way, they offer reliable framework for study of the notorious 
lidar ratio problem (i.e., extinction to backscatter codependency), widely discussed 
in the pure lidar community [16]. Since, for example in the case of the Pierre Auger 
Observatory, where optical properties have to be known over large volumes of at- 
mosphere, and a scanning lidar is therefore a necessity, both mentioned methods 
represent natural first choice of data analysis. 
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